arXiv: 1508.04781 v 1 [astro-ph.HE] 19 Aug 2015 


PROCEEDINGS 

0F SCIENCE 



Observation of Anisotropy in the Arrival Direction 
Distribution of TeV Cosmic Rays with HAWC 


Segev Y. BenZvi, 0 Daniel W. Fiorino^ and Stefan Westerhoff' for the HAWC 
Collaboration 0 

a Department of Physics and Astronomy, University of Rochester, Rochester, NY, USA 
b Department of Physics and WIPAC, University ofWisconsin-Madison, Madison, WI, USA 
c For a complete author list, see www.hawc-observatory.org/collaboration/icrc2015.php. 

E-mail: fiorino@wipac.wisc.edu 

The High-Altitude Water Cherenkov (HAWC) Observatory, located 4100 m above sea level near 
Sierra Negra (19° N) in Mexico, is sensitive to gamma rays and cosmic rays at TeV energies. 
The arrival direction distribution of cosmic rays at these energies shows significant anisotropy on 
several angular scales, with a relative intensity ranging between 10 -3 and 10 -4 . We present the 
results of a study of cosmic-ray anisotropy based on more than 86 billion cosmic-ray air showers 
recorded with HAWC since June 2013. The HAWC cosmic-ray sky map, which has a median 
energy of 2 TeV, exhibits several regions of significantly enhanced cosmic-ray flux. We present 
the energy dependence of the anisotropy and the cosmic-ray spectrum in the regions of significant 
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1. Introduction 

The HAWC Observatory has recently completed construction at 4100 m above sea level near 
Puebla, Mexico. The observatory, located at 19°N latitude, is designed to study the sky in gamma 
rays and cosmic rays between 50 GeV and 100 TeV. While cosmic rays are the major source of 
background in the gamma-ray analysis, the distribution of the arrival directions of the cosmic rays 
is itself of significant interest. Cosmic rays with rigidities of several TV strongly scatter on Galactic 
magnetic fields (r g « 0.001 pc), scrambling their initial direction many times before reaching Earth. 
Despite this, the last decade of experiments has reported stable 10 -3 — 10 -4 sidereal anisotropy in 
the arrival direction distribution of TeV cosmic rays in both hemispheres (see Ref [1] and references 
therein). The anisotropy has been observed at large (> 60°) and small angular scales by multiple 
experiments. 

At TeV energies, the large-scale structure is dominated by dipole and quadrupole moments 
of consistent phase throughout all experiments. These studies require year-long observations to 
cancel the effects of the Solar dipole. This study will not be presented with this half-year data set. 
After removal of the large-scale structure, significant small-scale structure appears, most notably a 
region of excess near (a = 60°, 8 = —5°) [2, 3]. Small-scale deficits do not appear to be significant 
features and may be an artifact of the removal of the large-scale structure. 

The origin of the anisotropy is not well understood. It has been suggested that weak dipole or 
dipole-like features should be a consequence of the diffusion of cosmic rays from nearby sources in 
the galaxy [4, 5]. This hypothesis is currently of considerable interest [6,7, 8]. Excesses may be the 
result of magnetically connected cosmic-ray accelerators or something more exotic [9, 10]. Small- 
scale anisotropy may also arise as a feature of turbulent scattering of large-scale anisotropy [11, 12]. 

The newly finished HAWC detector observes a larger portion of the sky with a higher recorded 
event rate than any high-energy cosmic-ray detector before it. This proceeding serves as an update 
to Ref. [1]. These data cover a region of the sky previously unobserved by experiments operat¬ 
ing in the northern and southern hemispheres. In these proceedings we present measurements of 
anisotropy on small angular scales with an emphasis on the energy-dependence of the significant 
features. We compare the observed anisotropy with previous measurements of the northern and 
southern skies. 

2. Data Set 

The HAWC Observatory is a 22,000 m 2 array of close-packed water Cherenkov detectors 
(WCDs). Each WCD consists of a cylindrical steel water tank 4.5 m in height and 7.3 m in diameter. 
The complete detector now comprises 300 WCDs with 1200 PMTs. 

The analysis presented in this paper uses the data collected during the operation of 111 WCDs 
(HAWC-111) between June 16, 2013 and July 8, 2014. Triggers for gamma-ray and cosmic- 
ray air showers were formed with a simple multiplicity trigger which requires >15 PMTs to be 
above threshold within a sliding time window of 150 ns. The trigger rate in HAWC-111 was 
approximately 15 kHz. 

A subset of the HAWC-111 data set was used for the anisotropy search. A cut in the fractional 
number of PMTs > 6% is used to remove poorly reconstructed events from the data. Additionally, 
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only continuous sidereal days of data were chosen for in order to reduce the bias of uneven exposure 
along right ascension. The resultant data set has a livetime of 181 days with 8.6 x 10 10 well- 
reconstructed cosmic-ray air showers. Using the detector simulation we estimate that the median 
energy of the data set is about 2 TeV. We estimate the angular resolution for the primary particle 
direction to be approximately 1.0° from both simulation and observation of the cosmic-ray Moon 
shadow. This is sufficient to observe the 5° — 180° features in the anisotropy of the cosmic rays. 

3. Analysis 

To search for anisotropy, we directly compute the relative intensity as a function of equatorial 
coordinates (a,<5). We begin by binning the sky into an equal-area grid with a resolution of 0.2° per 
bin using the HEALPix library. A binned data map N(a, 8)i is used to store the arrival directions 
of air showers recorded by the detector for each angular bin i. 

To produce residual maps of the anisotropy of the arrival directions of the cosmic rays, we 
must have a description of the arrival direction distribution if the data arrived isotropically at Earth, 
(N) (a, 8)i. We calculate this expected flux from the data themselves in order to account for all-sky 
rate variations and the changing viewing angle of each pixel. The spatial distribution of events is 
unique to a detector configuration and stable for periods exceeding 24 hours. The all-sky rate varies 
in Solar time by ~ 5% and approximates a sinusoid. 

The reference map (A) (a, 5)/ is produced using the direct integration technique described in 
Ref. [13], adapted for the HEALPix grid. In brief, we proceed by collecting all events recorded 
during a predefined time period At and integrate the local arrival direction distribution against the 
detector event rate. The method effectively smooths out the true arrival direction distribution in 
right ascension on angular scales of roughly At • 15° hr -1 such that the analysis is only sensitive 
to structures smaller than this characteristic angular scale. The direct integration procedure also 
compensates for variations in the detector rate while preserving the event distribution in declination. 
Once the reference map is obtained, we calculate the deviations from isotropy by computing the 
relative intensity 
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which gives the amplitude of deviations from the isotropic expectation in each angular bin i. The 
significance of the deviation can be calculated using the method of Li and Ma [14]. 

This analysis method can be sensitive to any maximum angular scale through the choice of 
At. Due to the sampling of the reference map along lines of right ascension, the maximum angle 
scales as (cos<5) -1 . Only a choice of 24 hours ensures a uniform angular scale as a function of 
declination. To remove larger structure, a multipole fit can be subtracted to access lower angular 
scales while preserving the maximum angular scale throughout the map. 


4. Measurement of the Small-scale Anisotropy 

The analysis was carried out on HAWC-111 data using A t- 24 hr to obtain sensitivity to all 
angular features equally over the sky. The results are plotted in Figure 1 in relative intensity. The 
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Figure 1: Relative intensity of the cosmic-ray flux for 181 days of HAWC-111, in equatorial coordinates. 
A period At = 24 hr is used to access the largest features present in the map. The map is shown with 10° 
smoothing applied. 


data have been smoothed using a 10° top-hat function to make the clustering of arrival directions 
readily apparent. At this time we do not provide an in-depth description of the observed large-scale 
features. The measurement in Fig. 1 is a combination of sidereal anisotropy and the Solar dipole 
effect which causes an excess of cosmic rays in the direction of the Earth’s motion around the 
Sun. The Solar contamination in this data set has been confirmed by observing significant signal 
in the unphysical coordinate system using “anti-sidereal" time. The two are difficult to disentangle 
without a full year of data in which the signals make a complete transit in the other reference frame. 

The large-scale signals can be subtracted from the map in Fig. 1 using a multipole fit to all 
multipole moments £ < £ max . Any small-scale residual from the Solar signal should be negligible 
(< 10 -5 ), so we are left with the small-scale structure in the sidereal frame. The residual structure 
will contain power at angular scales less than the 180°/^ m<3X . To reveal angular scale < 60°, a 
choice of £ max = 3 was made. The fit of all moments with £ < 3 in the map in Fig. 1 were was 
subtracted to produce the relative intensity and significance maps shown in Figure 2. Again, a 10° 
smoothing has been applied, so only pre-trials significances are shown. The estimated number of 
trials is at most equal to the number of independent bins 10 5 ), but is in fact much smaller, since 
we are not performing a blind search; all excess regions have been perviously observed by other 
experiments. The subtraction of the fit removes angular power in £ < 3 without influencing the 
higher order multipoles as evidenced by a comparison of the angular power spectra in Figure 3. 

With the trials factor taken into account, three significant features remain, the strong region 
of excess flux at (a = 60°, 8 = —10°), the band of excess along a = 120°, and a weaker excess 
at (a = 240°, 8 = 45°). These excesses lie coincident with the 10° — 20° regions of cosmic-ray 
access first observed by Milagro (Regions A and B in [2]) and a third region observed only by 
ARGO-YBJ [3] (Region C). The strongest excess, with a pre-trial significance of 20.6c, is found at 
a = 58.4° and 8 = —8.8° and corresponds to Region A in the Milagro map. The relative intensity 
of the excess in this region peaks at (7.9 ± 0.4 ± 0.4) x 10 -4 , where the first error is statistical 
and the second error is systematic. The median cosmic-ray energy at this declination is 2.1 TeV. 
For comparison, we fit a 2-D Gaussian function to the relative intensity map around Region A. The 
center is located at a = 59.2° ±0.6° and 8 = —7.2° ±0.7°, with an amplitude of (8.7 ±0.9) x 10 -4 . 
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Figure 2: Relative intensity (top) and pre-trial significance ( bottom ) of the cosmic-ray flux after fit and 
subtraction of the dipole, quadrupole, and octupole terms from the map shown in Fig. 1. The map is shown 
with 10° smoothing applied. 



Figure 3: Angular power spectra of the unsmoothed relative intensity map before (blue) and after (red) 
subtracting the large-scale structure (£ < 3). Gray errors bands show the 68% and 95% spread of the Q 
for isotropic data sets. Comparing the band to the data shows which el /-modes significantly contribute to 
the sky map. The error bars on the Q are the square root of the variances returned by a fit using a power 
spectrum estimator (PolSpice). 
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The width is nearly symmetric in a and <5, approximately 7.6° ± 1.2°. 

The elongated excess around a — 120°, coincident with Region B, extends over a wide range 
of declinations. It is most significant at (119.5°, 38.9°) with a pre-trial significance of 15.6c and 
a relative intensity of (5.3 ±0.4 ±0.2) x 10 -4 . A third excess region, Region C, is centered at 
a = 206.0° and 8 = 23.8° with a pre-trial significance of 10.5c and a peak relative intensity of 
(2.8 ± 0.3 ± 0.7) x 10 -4 . The median cosmic-ray energy at this declination is 2.0 TeV. 

There are also regions of strong deficits visible, typically on both sides of the strong excess 
regions. The appearance of these deficit regions, correlated with the excess regions, is likely an 
artifact caused by imperfections in the multipole fit. 

5. Study of Region A 

An early version of energy estimation has been developed for studying the HAWC data set as 
a function of energy. Seven energy bins are separated by cutting on both the incident zenith angle 
and the fraction of the detector that observed a cosmic-ray air shower. The energies of these bins 
are determined from simulation. At this time, the energy bins are heavily overlapping, still, a clear 
evolution of Region A as a function of energy is identifiable. This region exhibits a harder spectrum 
than background. 

We estimate the statistical significance of the hardness the in Region A spectrum by comparing 
the slope of a linear fit of 81 versus log (2?) in Figure 4 to similar fits performed at many random 
locations in the field of view. The random locations excluded a 15° circle centered on Regions 
A, B, and C. The distribution of slopes for the random locations follows a Gaussian distribution 
centered at zero with a width of 0.9 x 10 -4 . The best fit slope at the position of Region A is (3.8 ± 
1.1) x 10 -4 , 4.2<t away from the mean. The reduced X 2 of a fit to a sloped-line is 1.16 compared 
to 5.66 for a horizontal line. This further supports a spectrum that is harder than background. 

A zoomed view of Region A in four energy-proxy bins (the final four bins of Fig. 4 were com¬ 
bined) is provided in Figure 5. In addition to the spectral steepening, this figure shows the excess 
extending to higher declination as a function of energy, consistent with the ARGO observation. 

6. Discussion 

This half-year HAWC-111 data set is already one third as large as the Milagro 8-year data set 
and one fourth of the ARGO-YBJ 5-year data set. We have found significant small-scale structure 
that coincides with 3 of the 4 previously published regions of cosmic-ray excess. In general, there 
is excellent agreement between the observed excesses and deficits of previous measurements by 
ARGO-YBJ and Milagro in the Northern Hemisphere. 

Region A at a = 60° appears to have two sub-regions in the TeV range. At higher energies 
(>10 TeV), the northern part at 8 = 5° appears. Below this, the southern part at 8 = —5° is dom¬ 
inant. Milagro saw the northern part of Region A, but only observed down to 8 = —5°. There 
is a slight indication of the southern part of Region A in the unpublished relative intensity maps. 
Both the HAWC and ARGO maps see the change with energy in Region A. This suggests that the 
hardening of the spectrum of Region A is coincident with a northernly movement of the excess or 
the emergence of an entirely new region of excess. 
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Figure 4: Spectrum of Region A in relative intensity in different energy-proxy bins. The energies of the 
data were determined from simulation. The error bars on the median energy values correspond to a 68% 
containing interval. 
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Figure 5: Relative intensity of Region A for 4 different energy-proxy bins. The square mark denotes the 
location of the centroid of Region A as reported by Milagro (a = 69.4°, <5 = 13.8°) at 10 TeV. The median 
energy of the data in each plot from left to right is 1.7^3 TeV, 3.2+\° 4 9 TeV, TeV, and 14.TeV. 


Region B at a = 120° is visible at all declinations observed by HAWC and seems to follow 
the edge of the large-scale sidereal deficit. The extension of the excess to higher declinations is 
notable because this is not observed to high significance in the Milagro data, even though Milagro 
was located at a higher latitude (35°N). The significance maps published by ARGO also extend to 
high declinations, and these maps indicate that the Region B excess extends northward. 

Region C at a = 205° is not significant in the Milagro data, but the ARGO collaboration has 
observed a hot spot at the same location of similar shape and intensity. ARGO also observed a 
fourth significant region at a = 280° which appears as a sub threshold excess in the HAWC map. 
Continued observation with HAWC will help confirm this as a true excess. 

A comparison between data from HAWC-111 and IceCube (20 TeV median) in the Southern 
Hemisphere is very interesting because connecting the northern and southern measurements would 
eliminate biases from partial sky coverage. No clear connection of the small-scale anisotropy that is 
present in both hemispheres has been made yet. The published IceCube maps are of higher energy 
than HAWC and other northern measurements. Because of this, a first effort at combining HAWC 
and IceCube data uses cuts to bring their median energies closer to a central value [15]. 
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Ten years after the first sky maps that charted the cosmic-ray anisotropy, the state of the art 
in this field has moved towards studying significant cosmic-ray features as a function of rigidity. 
Precision measurements of cosmic-ray anisotropy are readily obtained from air shower detectors, 
such as HAWC. Initial improvements to the resolution and range of measured cosmic-ray energies 
using these detectors have revealed energy-dependent behavior of the anisotropy. Further studies 
of this phenomenon may provide insight into local cosmic-ray production and propagation or may 
even become a tool for studying exotic physics. 
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